Analysis of power simulations: actual data

Purpose of the document

In this document, we analyze the results of our simulations.

Loads

First, we load the results of our simulations and useful packages.

Quick analysis

First, we graph of the evolution of type M, type S and power for each quasi-experiment as a function of different parameters, holding other parameters constant. This function takes as input a “summary_simulations” type of data frame.

graph_evol_by_exp <- function(df, var_param = "n_days_study", stat = "power") {
  
  var_param_name <- str_replace_all(var_param, "_", " ")
  stat_name <- str_replace_all(stat, "_", " ")
  
  #considering baseline values
  df_filtered <- df %>% 
    filter(str_detect(formula, "temperature")) #to only consider the model with all covariates
  
  if (var_param != "p_treat") {
    df_filtered <- df_filtered %>% 
      filter(p_treat == sim_param_base[["p_treat"]])
  } 
  if (!(var_param %in% c("n_days_study", "average_n_obs"))) {
    df_filtered <- df_filtered %>% 
      filter(n_days_study == sim_param_base[["n_days_study"]])
  } 
  if (var_param != "percent_effect_size") {
    df_filtered <- df_filtered %>% 
      filter(percent_effect_size == sim_param_base[["percent_effect_size"]])
  }
  
  #graph itself
  graph <- df_filtered %>% 
    mutate(
      quasi_exp = str_to_sentence(str_replace_all(quasi_exp, "_", " "))
    ) %>% 
    ggplot(aes(x = .data[[var_param]], y = .data[[stat]])) + #, color = .data[[quasi_exp]] + 
    geom_point() +
    geom_line(linetype = "dashed", size = 0.1) +
    facet_wrap(~ quasi_exp) +
    ylim(c(0, ifelse(stat == "power", 100, NA))) +
    labs(
      title = paste(
        str_to_title(stat_name), ifelse(stat == "power", "increases", "decreases"),
        "with", var_param_name
      ),
      subtitle = "Comparison across quasi-experiments",
      x = var_param_name,
      y = str_to_title(stat_name)
    ) 
    # theme(legend.position = "none")
  
  return(graph)
} 

# graph_evol_by_exp(summary_simulations)

We then plot all the graphs. To do so, we create a tibble containing every parameter and statistics we want to plot and map our function on this tibble.

stat <- c("power", "type_m", "type_s")
var_param <-  c("n_days_study", "average_n_obs", "percent_effect_size", "p_treat")

param_by_exp <- crossing(stat, var_param)

pmap(param_by_exp, graph_evol_by_exp, df = summary_simulations)

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

We can also plot the graphs for each parameter separately, for a better readability.

Length study

map(
  stat, 
  graph_evol_by_exp, 
  df = summary_simulations, 
  var_param = "n_days_study"
)

[[1]]

[[2]]

[[3]]

map(
  stat, 
  graph_evol_by_exp, 
  df = summary_simulations, 
  var_param = "average_n_obs"
)

[[1]]

[[2]]

[[3]]

Effect size

map(
  stat, 
  graph_evol_by_exp, 
  df = summary_simulations, 
  var_param = "percent_effect_size"
)

[[1]]

[[2]]

[[3]]

Proportion of treated days

map(
  stat, 
  graph_evol_by_exp, 
  df = summary_simulations, 
  var_param = "p_treat"
)

[[1]]

[[2]]

[[3]]